AD-A268  768 


Office  of  Naval  Research 


Grant  N00014-90-J-1871 


Technical  Report  No.  14 

THE  LONGITUDINAL  FREE  VIBRATIONS 
OF  A  FIXED-FREE  TWO-PHASE  ELASTIC  BAR 


by 

Jim  M.  Restuccio 


Division  of  Engineering  and  Applied  Science 
California  Institute  of  Technology 
Pasadena,  California  91125 


Thif  iorunirni  has  Dean  approved 
for  putli-  i.nd  sale;  its 

d;-tr;r  -tiO;:  uriiinited. 


June,  1993 


93**  19997 


8  26  00  8 


1 


1.  Introduction 


In  this  paper,  the  longitudinal  free  vibrations  of  a  fixed-finee  bar  are  studied.  It 
is  assumed  that  the  bar  initially  consists  of  two  phases,  one  of  which  was  obtained 
from  the  other  by  a  maitensidc  phase  transformation.^  It  is  also  that  both 

phases  of  the  bar  have  elastic  constitutive  behavior.  For  a  bar  that  consists  entirely 
of  one  phase  that  behaves  elastically,  it  is  well  known  that  during  the  fiee  vibnUions 
of  the  bar  the  displacement  and  stress  at  each  point  of  the  bar  oscillate  as  time 
progresses  [4].  If  there  is  damping  present,  these  oscillations  will  decay  and  go  to 
zero  as  dme  goes  to  infinity,  otherwise  their  amplitudes  will  remain  constant  in  dme. 

G)nsidering  this,  for  a  bar  that  inidally  consists  of  two  different  phases  that  both 
behave  elasdcally,  one  might  expect  that  the  displacement  and  stress  at  each  point 
of  the  bar  will  also  have  oscillatory-type  behavior  during  the  free  vibradons  of  the 
bar.  If  this  is  the  case,  the  driving  tracdon  at  the  interface  separating  the  two  phases 
will  oscillate.  As  a  result  of  this,  if  the  nominal  phase  boundary  velocity  is  related 
to  the  driving  tracdon  through  a  kinedc  rel^on  that  does  not  have  an  interval  of  the 
driving  tracdon  corresponding  to  a  zero  nominal  phase  boundary  velocity,  the  nomiiud 
phase  boundary  velocity  will  also  oscillate.  Since  energy  is  dissipated  when  the  phase 
boundary  moves  and  passes  over  pardcles  of  material  of  one  phase  converting  them 
into  pardcles  of  material  of  the  other  phase,  one  might  conclude  that  the  oscillatory- 
type  response  of  the  two-phase  bar  during  the  free  vibradons  of  the  bar  should  decay 
as  dme  increases.  It  is  this  damping  behavior  of  the  two-phase  bar  that  will  be  the 
main  subject  of  this  chapter.  The  solutions  of  the  boundary  value  problem  will  be 
determined  by  a  numerical  method,  and  the  damping  of  the  two-phase  bar  will  be  ’  ^ 

■  —  ■  — ■  .  j 

'  ManensiticphaieiiaiisfonDatioiu  ate diffusioolesi<olid-i(dkipbawtnBsfoaDatiafiiwfaidi,>iaaiiso(lier things,, 

have  condnuoua  displacements,  with  possible  discootinuons  siiains,  at  thek  phase  boundary.  These  types  rf  -  —  :  ~ . . . 

tiansfonnatioiis  ate  also  chaiacteriied  by  the  product  phase  bavins  a  shape  defonnatkmielative  to  the  nadefoaned  doL 

parent  phase,  which  cotrespoods  to  an  unstressed  undefonned  oonfigunlion  of  that  product  phase.  See  [S]  andV^ -r.  o  f  v. 
[8]  for  more  detailed  and  comprehensive  discussions  about  these  types  of  tansfixmadoos.  ' 
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Studied  as  the  material  coefficients  are  varied.  The  values  of  the  material  coefficients 
that  result  in  the  maximum  damping  will  also  be  investigated. 

2.  The  Kinematics 

A  one-dimensional  finite  bar  that  initially  consists  of  two  phases  is  considered.  It 
is  assumed  that  the  process  under  consideration  occurs  in  a  time  interval  F  =  [toi  fi]. 
Additionally,  for  the  problem  that  is  considered,  a  continuum  model  that  was 
developed  in  [6]  (see  also  [7])  is  used. 

G>nsider  a  stationary  reference  configuration  R  for  the  bar.  Let  x  denote  a  point 
in  R  and  let  Z.  be  the  length  of  the  bar  with  respect  to  this  reference  configuration. 
Considering  this,  R  can  be  expressed  as  R  =  {xf  x  e  [0,X]}.  Let  y(x,t)  be  the 
suitably  smooth  and  invertible  mapping  which  maps  R  into  the  deformed  configuration 
of  the  bar  at  each  t  €  F,  with  y(x,t)  —  x  +  u(x,<)  V  (x,t)  €  R  x  F.  The  quantity 
u(x,  t)  represents  the  displacement  of  a  particle  of  material  at  y  =  y(x,  t)  from  the 
point  X  6  R  at  time  t  €  F.  In  the  following,  the  two  phases  of  the  bar  will  be 
referred  to  as  phase  1  and  phase  2.  Let  x  =  s{t)  be  the  reference  position  of  the 
phase  boundary  separating  phase  1  from  phase  2  at  time  t  €  F.  It  is  assumed  that 
particles  of  material  with  reference  points  in  R~  =  {x/  x  €  [0,s(<)]}  at  time  t  €  F 
are  in  phase  1,  and  it  is  assumed  that  particles  of  material  with  reference  points 
in  R*^  =  {x/  X  €  [s(t),I]}  at  time  t  €  F  are  in  phase  2.  It  is  assumed  that  R~ 
coincides  with  an  unstressed  undeformed  configuration  of  phase  1.  We  next  assume 
that  there  exists  a  shape  deformation  of  phase  2  with  respect  to  K*"  that  corresptHids 
to  an  unstressed  undeformed  configuration  of  that  phase.  Let  R^  be  the  referraoe 
configuration  coinciding  with  this  shape  deformation  for  all  t  €  F.  Let  xi  denote  a 
point  in  Rf,  and  let  Xi(x,t)  be  the  suitably  smooth  and  invertible  mapping  which 
maps  R'*'  into  R?"  at  each  t  €  F,  with  X|  =  xi(x,  t)  V  x  €  R"*"  at  each  t  G  F.  For  the 
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following  problem,  it  is  assumed  that  £i  is  given  by 

x,(x,f)  =  x+7o(a:-s(t)).  (2.1) 

The  displacement  gradient  (transformation  strain)  corresponding  to  this  choice  of  X| 
is  70,  the  Jacobian  is  J  =  1  +  70,  and  Rf  =  {xi/xi  e  [xi(s(<),f),Xi(L,f)]}  = 
{xi/xi  €  [s(f),Xi(f)]},  where  L\{t)  =  xi(X,f)  =  L  +  7o(X-s(f))-  It  is  also 
assumed  that  70  >  —1  so  that  reflections  are  excluded  from  (2.1).  Additionally,  flie 
inverse  z  of  the  xt  given  by  (2.1)  is 

x(z,,  t)  =  — (zi  +  7os(<)),  (2.2) 

1+70 

V  xi  €  at  each  t  e  F.  The  mapping  which  maps  Rf  into  the  deformed 
configuration  of  phase  2  at  each  t  6  F  is  represented  by  yi(xi,f),  with  yi(zi,t)  s 
xi  +  ui(zi,<)  V  xi  €  Rf  at  each  <  €  F. 

3.  The  Continuum  Model 

In  the  continuum  model  that  is  used  here,  the  constitutive  equations  for  each 
phase  are  defined  with  respect  to  different  reference  configurations.  More  specifically, 
the  constitutive  equations  for  phase  1  are  defined  with  respect  to  R~,  and  the 
constitutive  equations  for  phase  2  are  defined  with  respect  to  R^.  Additionally,  the 
field  equations  for  phase  1  are  expressed  with  respect  to  R~,  and  the  field  equations 
for  phase  2  are  expressed  with  respect  to  R^  (see  [6],  [7]).  The  main  advantage 
of  using  this  continuum  model  for  the  problem  under  consideration  is  that  the  field 
equations  are  in  forms  that  permit  direct  linearization.  This  is  the  case  since  the 
displacements  for  each  phase  are  measured  from  a  reference  configuration  coinciding 
with  an  unstressed  undefoimed  configuration  of  that  phase,  and  consequently,  for 


4 


the  appropriate  boundary  and  initial  conditions,  the  displacement  gradients  can  be 
considered  intinitesimal.^ 


4.  The  Field  Equations  and  Jump  Conditions 


It  is  assumed  that  the  process  under  consideration  is  a  purely  mechanical  process 
with  no  body  forces  present  The  general  field  equations  and  junq>  conditions  using 
the  type  of  continuum  model  described  in  the  previous  section  and  for  a  purely 
mechanical  process  were  derived  and  discussed  in  [6],  [7]. 

The  field  equations  for  the  problem  under  consideration  consist  of  the  balance 
of  linear  momentum  for  phase  1  and  the  balance  of  linear  momentum  for  phase  2. 
These  equations  are 


da 


(4,1) 


V  *  €  R"  at  each  t  €  F,  and 


dci 

dxi 


Pi^i, 


(4^) 


V  xi  €  Rf  at  each  t  €  F,  respectively,  where  p{x)  is  the  density  of  the  material  per 
unit  volume  of  R,  pt(X|,t)  s  p{x)fJ  is  the  density  of  phase  2  per  unit  volume 
of  Rf,  is  the  nominal  stress  with  respect  to  R~  for  phase  1,  (Ti(z|,t) 

is  the  nominal  stress  with  respect  to  Rf  for  phase  2,  a(z,t)  =  ^y(z,f),  and 

As  mentioned  previously,  the  displacements  at  a  phase  boundary  separating  two 
phases  involved  in  a  martensitic  phase  transformation  are  continuous,  while  the  strains 
may  be  discontinuous.  Considering  this,  we  require  that  y  be  continuous  and  the  first 


^  It  is  also  sssumod  here  (hat  the  unstressed  undefonned  coofiguiatkw  of  each  phase  conespoods  to  a  relative 
minanum  of  (he  elastic  potential  for  that  phase. 
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and  second  derivatives  of  y  be  piecewise  continuous  on  R  x  F,  with  discontinuities 
occuiing  only  at  x  =  s(t).  The  continuity  of  displacement  condition  in  its  most 
direct  form  is 


y+(s(f),t)  =  y-(s(t),t).  (43) 

For  the  Xi  given  by  (2.1),  (4.3)  reduces  to 

u^(s(t),t)^u-(s(t),t).  (4.4) 

Considering  the  continuity  requirements  on  y  and  anticipating  the  form  of  the 
constitutive  equations  for  the  type  of  material  under  consideration,  we  require  that 
the  stress  be  piecewise  continuous,  with  discontinuities  occuring  only  at  x  =  s(f). 
The  jump  condition  at  x  =  )  representing  the  balance  of  linear  momentum  is 

p(ij^  -  v~)i  =  0,  (43) 

where  i(<)  =  Differentiating  (4.3)  with  respect  to  time  yields 

(to  +  7^70  +  7^  -  7")s  +  =  0,  (4.6) 

at  X  =  a(f),  where  7  =  ^  and  7i  =  Using  (4.6)  in  (4.5),  we  can  obtain  an 
alternate  form  for  the  linear  momentum  jump  condition: 

=  p(yo  +  7^70  +  7^  "  7")(«)^-  (4.7) 

The  remaining  equation  at  the  phase  boundary  is  a  kinetic  relation  relating  i  and 
the  driving  traction  /  (see  [1]).  This  kinetic  relation  is  a  constitutive  equation  and 
will  discussed  in  the  next  section.  For  the  problem  under  consideration  the  driving 
traction  is  given  by 

/  =  J  i  (<7^  +  O’” )  (70  +  7^70  +  7^  -  7"  )  • 


(4.8) 
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We  also  have  the  boundary  conditions  at  x  =  0  and  xi  =  Li(t)  (or  x  =  L).  For 
the  fixed  boundary  condition  at  x  =  0,  we  have 

u(0,0  =  0,  (4.9) 

Vt  €  r.  At  the  xi  =  Li(i)  boundary,  there  is  nothing  applied  to  th^  boundary;  i.e.  it 
is  a  £tee  boundary.  Therefore,  the  traction  at  this  boundary  is  necessarily  zero.  Thus, 
the  boundary  condition  at  xi  =  Li(i)  is 

at(X,(i),t)  =  0,  (4.10) 


Vt  €  r. 


5.  The  Constitutive  Equations 

It  is  assumed  that  both  phase  1  and  phase  2  are  homogeneous  elastic  (or 
hyperelastic)  phases.  In  particular,  for  phase  1  we  assume  that  there  exists  an  elastic 
potential 

W  =  Wiy),  (5.1) 


such  that 


cr  = 


dW 
d7  ’ 


(52) 


and  for  phase  2  we  assume  that  there  exists  an  elastic  potential 


W^  =  ir,(7i), 


(53) 


such  that 


(5.4) 
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If  we  solve  for  ui  =  ui(xi,t)  in  the  dynamic  boundary  value  problem,  the 
acceleration  term  di  in  Equation  (4.2)  will  contain  several  ineitial-type  terms  that 
are  solely  a  result  of  Xi  being  a  function  of  time.  This  was  discussed  for  the 
general  three-dimensional  problem  in  [6],  [7].  As  was  further  discussed  there,  but 
again  for  the  general  three-dimensional  problem,  these  ineitial-type  terms  can  be 
avoided  by  instead  solving  for  U|  =  tii(x,t)  in  the  boundary  value  problem,  where 
ui(x,<)  =  ui(xi(x,t),t)  and  ui(xi,<)  =  tii(x(xi,f),t).  This  will  be  done  in  the 
boundary  value  problem  that  is  considered  here. 

It  is  assumed  that  phase  1  is  unstressed  at  7  =  0  and  phase  2  is  unstressed  at  71  = 
0.  We  next  assume  that  the  initial  conditions  are  such  that  I7I  «  1  V  x  €  [0,  s(t)) 
and  |7i|  «  1  V  X  €  (s(t),I],  at  each  t  €  F,  where  here  and  in  the  following, 
7,  =  a:  ^  For  the  x  given  by  (2.2),  71  »  For  these  assumpdons, 

W  for  phase  1  can  be  written  as 


(5.8) 
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We  note  that  the  constitutive  quantities  for  phase  1  given  by  (5.1),  (5.2),  (5.5),  and 

(5.7)  are  defined  with  respect  to  R~,  and  the  constitutive  quantities  for  phase  2  given 
by  (5.3),  (5.4),  (5.6),  and  (5.8)  are  defined  with  respect  to  Rf,  even  though  (5.6)  and 

(5.8)  are  in  terras  of  the  independent  variable  x  e  R'*'. 

We  next  assume  that  W’  =  J  W".  This  assumption  might  be  most  ^propriate 
if  phase  1  and  phase  2  represent  two  different  variants  of  the  same  martensite.  We 
note,  however,  that  if  this  were  the  case  we  need  not  assume  that  Ei  =  E,  since 
in  a  real  three-dimensional  material  the  moduli  in  a  given  direction  of  two  variants 
of  martensite  separated  by  a  phase  boundary  are  not,  in  general,  the  same,  and  this 
can  be  incorporated  into  a  one-dimensional  model  by  assuming  that  Ei  ^  E.  For 
these  assumptions  and  using  (5.5)-(5.8),  the  first-order  ^proximation  of  the  driving 
traction  given  by  (4.8)  is 


/« 


(5.9) 


We  next  postulate  a  kinetic  (constitutive)  relation 

i  =  W),  (5.10) 

such  that 

$(/)/  >  0,  (5.11) 

for  all  /  [1].  The  requirement  given  by  (5.1 1)  is  imposed  so  that  energy  is  dissipated 
(or  preserved  if  the  equality  sign  holds)  during  a  phase  transformation,  instead  of 
being  created.  We  further  assume  that  for  the  problem  under  consideration,  the  first- 
order  form  of  $(/)  is 
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(5.12) 


where  the  constant  v  depends  on  the  material  and  is  such  that  1/  >  0  so  that  (S.ll) 
is  satisfied. 


6.  The  Boundary  Value  Problem 

Substituting  (5.7)  into  (4.1)  and  expressing  the  acceleration  in  that  equation  in 
terms  of  u,  we  obtain  the  following  for  the  balance  of  linear  momentum  for  phase  1: 


(6.1) 


Va:  €  (0,3(t))  at  each  t  6  F.  Substituting  (5.8)  and  pi  =  p/ J  =s  p/(l  +70)  into 
(4.2),  and  calculating  di  in  that  equation  using  yi(xi(x,t),t)  =  X|(x,t)  +  ui(x,t), 
we  obtain  the  following  for  the  balance  of  linear  momentum  for  phase  2: 


(I  +  70)  5x2  Qti  )' 


{62) 


V  X  G  (s(<),  L)  at  each  t  €  F, 

For  the  linearized  problem,  the  continuity  of  displacement  condition  is  still  given 
by  (4.4),  since  there  is  nothing  to  linearize  in  that  equation.  We  next  consider  the 
linearized  form  of  the  linear  momentum  juiiq>  condition  given  by  (4.7).  Because  i  is 
related  to  7  and  71  through  the  Idnefic  relation  given  by  (5.12),  is  second-tMrder  in 
7  and  7i ;  i.e  goes  to  zero  at  the  same  rate  as  7^,  7,^,  and  771  go  to  zero.  However, 
the  second-order  terms  in  also  contain  the  constants  and  EE\lv. 
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Considering  the  fact  that  the  magnitudes  of  E  and  E\  are  very  large,  if  the  magnitmig 
of  1/  is  not  also  large,  7  and  71  might  have  to  be  unrealistically  small  in  order  for  the 
terms  in  to  have  magnitudes  that  are  negligible  in  comparison  to  the  magnitudes 
of  the  first-order  terms.  In  the  following,  it  will  be  assumed  that  the  values  of 
the  material  coefficients  are  such  that  the  magnimdes  of  the  second-order  terms  in 
Equation  (4.7)  are  negligible  in  comparison  to  the  magnitudes  of  the  first-order  terms 
in  that  equation  for  realistic  values  of  the  infinitesimal  strains.^  The  restrictions  that 
this  assumption  puts  on  the  relative  values  of  u  and  the  other  material  coefficients  can 
best  be  observed  when  the  linear  momentum  jump  condition  and  the  kinetic  relation 
are  in  nondimensional  form,  which  will  1 done  in  the  next  section.  The  true  first- 
order  s^proximadon  of  the  linear  momentum  jump  condition  given  by  (4.7)  is 


El 

ditt 

E 

du 

l+7o 

.  . 

(•(*).*) 

(•aw 


=  0, 


(63) 


Vt  €  r,  which  is  equivalent  to  the  continuity  of  tractions  across  the  interface.  Using 
(6.3)  in  (S.12),  the  first-order  approximation  of  the  kinetic  relation  for  the  problem 
under  consideration  can  be  written  as 


a  =  — 


70-E 


dx 


(6.4) 


J  (•(*),«) 


Vt  €  r. 


For  the  linearized  problem,  the  fixed  boundary  condition  is  still  given  by  (4.9),  and 
using  (S.8)  in  (4.10),  the  first-order  ^proximation  of  the  free  boundary  condition  is 

3ui| 


dx 


=  0, 


(63) 


(i.«) 


Vt  €  r. 


’  We  note  that  the  values  of  the  stnini  at  the  interface  as  time  progresses  are  ptopottkmal  to  the  initial  conditinns 
that  ate  giveiL 
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7.  The  Nondimensional  Form  of  the  Boundary  Value  Problem 

For  the  problem  under  consideration,  we  define  the  following  nondimensional 


variables: 


1  =  -,  t=Ljt,  &\i)  =  — 2 — ’ 


u{x,i)  = 


u(xL,i/u}) 


,  ui{x,t)  = 


Ui{xL,ifui) 


(1  +  7o)£’ 


yay/Ep' 


where  w  =  y/W/{pU).*  Using  these  nondimensional  variables,  the  nondimensional 
form  of  the  balance  of  linear  momentum  for  phase  1  is 

df^u  d^u 

IP-W'  (7^) 


Vx  €  (0,  J(^)  at  each  f  €  f ,  where  f  =  [tofi]  =  and  the  nondimensional 

form  of  the  balance  of  linear  momentum  for  phase  2  is 


Vx  €  (5(<),l)  at  each  f  €  f . 

The  nondimensional  form  of  the  continuity  of  displacement  condition  is 


V  t  €  r,  and  the  nondimensional  form  of  the  first-order  form  of  the  balance  of  linear 
momentum  jump  condition  given  by  (6.3)  is 


*  In  the  following,  x  will  denote  the  nondimensional  independent  variable  defined  by  (7.1)| ,  and  not  the  inverse 
function  of  the  function  xi  which  maps  K*  into  at  each  t  C  F. 
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V  f  €  f .  Additionally,  the  nondimensional  form  of  the  kinetic  relation  given  by 
(6.4)  is 


Vt  €  f .  In  the  previous  section,  the  issue  concerning  when  it  is  appropriate  to  neglect 
the  second-order  terms  in  the  linear  momentum  jump  condition  was  discussed.  The 
nondimensional  form  of  the  lowest-order  term  that  was  neglected  in  that  equation 
is  7o(^)  .  From  this,  (7.6),  and  the  definitions  of  the  nondimensional  parameters, 
we  can  observe  what  relative  values  of  the  material  coefficients  are  appropriate  for 
the  assumption  that  the  magnitude  of  7o(^/  is  negligible  in  comparison  to  the 
magnitudes  of  the  terms  in  (7.5),  for  realistic  values  of  the  initial  conditions. 

The  nondimensional  form  of  the  fixed  boundary  condition  is 

u(0,t)  =  0,  (7.7) 

V  f  €  f ,  and  the  nondimensional  form  of  the  free  boundary  condition  is 


du,| 


=  0, 


(7.8) 


V  t  €  f . 
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We  must  also  specify  initial  conditions  for  the  two-phase  bar.  In  particular,  for 
phase  1  we  specify 


u(z,<o)  =  A(x), 


du 


(7.9) 


for  0  <  X  <  s(t),  and  for  phase  2  we  specify 


tit  (x,  to)  =  ^i(x), 


dui 

W 


=  ffi(x) , 


(7.10) 


for  i(t)  <  X  <  1.  We  also  must  specify  an  initial  position  for  the  phase  boundary. 
In  particular,  we  specify 


5  (to)  =  So- 


(7.11) 


The  initial  boundary  value  problem  that  will  be  considered  consists  of  the  field 
equations  (7.2)  and  (7.3),  the  continuity  of  displacement  condition  (7.4),  the  linear 
momentum  jump  condition  given  by  (7.5),  the  kinetic  relation  given  by  (7.6),  the  fixed 
boundary  condition  (7.7),  the  free  boundary  condition  (7.8),  and  the  initial  conditions 
given  by  (7.9H7.11). 

8.  The  Numerical  Method  of  Solution 

We  can  observe  that  the  differential  equations  involving  time  derivatives  in  the 
boundary  value  problem  presented  in  Section  7  consist  of  a  wave  equation  given  by 
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(7.2),  a  forced  wave  equation  given  by  (7.3),  and  an  ordinary  differential  equation 
given  by  (7.6).  The  boundaiy  conditions  for  Equations  (7.2)  and  (7.3),  given  by  (7.7) 
and  (7.8),  are  both  with  respect  to  fixed  boundaries.^  However,  the  jump  conditions 
given  by  (7.4)  and  (7.S),  which  are  also  types  of  boundary  conditions,  are  with  respect 
to  a  moving  boundary  (i.e.  the  interface).  We  also  note  that  because  s{i)  is  one  of  the 
unknown  dependent  variables  in  the  problem  and  the  jump  conditions  are  evaluated 
at  X  =  s(i),  the  boundary  value  problem  is  inherently  nonlinear  with  respect  to  s(t). 

Dynamic  boundary  value  problems  involving  moving  interfaces  within  finite 
bodies  have  been  studied  before  (see  [3]).  Some  of  the  more  well  known  of  these 
problems  are  the  class  of  problems  considered  to  be  Stephan  problems.  These  types  of 
problems  involve  melting  solids,  with  a  moving  interface  (or  boundaiy)  separating  the 
solid  from  the  liquid.  The  unknowns  in  these  types  of  problems  are  the  temperature 
distributions  of  both  the  liquid  and  solid  phases,  and  the  position  of  the  interface 
separating  these  two  phases.  The  governing  equations  for  these  Stephan  problems 
consist  of  heat  equations  for  both  phases  and  an  equation  governing  the  motion  of 
the  interface.  There  are  a  variety  of  numerical  methods  that  have  been  used  to  study 
these  Stephan  problems  (see  [3]  for  an  overview  and  discussion  of  dtese  methods). 
Among  these  numerical  methods  are  several  types  of  finite  difference  methods. 

For  the  problem  that  is  considered  here,  the  type  of  numerical  method  that  will  be 
used  is  a  finite  difference  method.  This  type  of  numerical  method  has  been  chosen, 
as  ^posed  to,  e.g.,  a  finite  element  method,  because  it  is  probably  the  most  straight 
forward  to  apply  to  the  type  of  boundary  value  problem  that  is  being  considered.  The 
particular  finite  difference  method  used  here,  however,  does  differ  somewhat  from 
the  finite  difference  methods  that  have  been  used  for  the  Stephan  problems  that  are 


’  Note  that  if  we  were  solving  for  «i  instead  of  «i  in  the  boundaiy  value  problem,  the  free  boundary  condition 
would  be  with  respect  to  a  moving  boundary. 
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discussed  in  [3].  Most  of  these  differences  reflect  the  fact  that  the  field  equations  in 
the  Stephan  problems  are  of  parabolic-type  with  monotonically  decaying  solutions, 
and  the  field  equations  in  the  problem  considered  here  are  of  hyperbolic-type  with 
decaying  oscillatory  solutions. 

In  the  following,  let  I[a,  h\  =  {m  ^  Zf  a  <m  <h,  ae  Z,  be  Z},  where  Z 
denotes  the  set  of  all  integers.  For  the  finite  difference  method  that  is  used  here, 
the  bar  will  be  divided  into  n  intervals,  each  of  length  h  =  1/n.  The  points 
I  =  th,  i  e  J(0,  n],  will  be  refeired  to  as  the  nodes  of  the  bar.  In  the  finite  difference 
method,  the  displacements  at  these  nodes  will  be  determined  (i.e  approximated).  The 
time  increment  is  denoted  by  T,  and  for  convenience  we  assume  in  the  following 
that  to  =  0.  Considering  this,  t  =  jT^  j  6  J[0,  /i],  where  ii  and  h  are  chosen 
such  that  liT  =  ti.  We  let  s(j)  represent  s(jT),  and  Jb(j)  denote  the  node  such 
that  j3(j)  -  A:(j)/i|  <  h/2,  at  time  i  =  jT.  Additionally,  we  let  p(j)h  represent  the 
distance  from  the  node  k{j)  -  1  to  a(jf),  at  time  i  —  jT.  Considering  this,  we  can 
write  s{j)  as  s(j)  =  [k(j)  —  1  at  each  j  e  l[0,/i].  Let  u(i,j)  represent 

u{ih,jT),  Q  <  ih  <  s(j)  at  each  j  €  1(0, /i],  and  let  Ui(»,i)  represent  u,(*A,jT), 
a(i)  <  ih  <  I  at  each  j  €  ^[0,/i].  Additionally,  let  q{j)  denote  the  displacement  at 
the  interface  at  time  t  =  jT;  i.e.  let  q{j)  represent  u(5(jT),jT)  =  u,(s(jT),jT). 
We  also  let  $(i)  represent  the  nondimensional  kinetic  relation  given  by  (7.6)  at  time 
i  =  jT. 

At  each  time  increment,  the  numerical  routine  begins  by  calculating  s{j  -f- 1). 
We  are  given  s(0)  as  an  initial  condition.  At  j  =  0,  we  use  an  Euler’s  method  to 
approximate  Equation  (7.6)  and  obtain  a(l).  In  particular,  we  use 

p(i) = m + |«(o) 


(8.1) 
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to  obtain  p(l).  Since  ib(0)  is  given  when  s(0)  is  given,  we  can  then  obtain  s(l).  Fot 
this  case  where  j  =  0,  we  can  use  the  initial  conditions  given  by  (7.9)i  and  (7.10)i  to 
obtain  the  derivative  term  in  $(0).  The  specific  form  of  these  initial  conditions  will 
be  discussed  at  a  later  point  in  this  section.  Also,  we  note  that  Euler’s  method  has 
an  0(T)  numerical  error.  At  j  =  1,  we  use  an  Adams-Bashfoith  two-step  method, 
which  has  an  O(T^)  error  (see  [2]).  to  approximate  Equation  (7.6)  and  obtain  s(2). 
The  resulting  equation  for  p(2)  is 

P(2)  =  P(l)  +  -  *(0)}  ,  (8.2) 

which  can  then  be  used  to  obtain  s(2).  Also,  to  obtain  the  derivative  term  in  $(1)  we 
can  use  the  initial  conditions  given  by  (7.9)  and  (7.10).  For  j  €  I[2,  li  -  1],  we  use 
an  Adams-Bashforth  three  step  method,  which  has  an  O(T^)  error,  to  approximate 
Equation  (7.6)  and  obtain  s(j  +  1).  In  particular,  we  use 

pU  +  1)  =  Pii)  +  -  16*0  -  1)  +  8*(j  -  2))  (8J) 

to  obtain  p(j  +  1)  for  j  €  I[2,/i  -  1].  p(j  +  1)  is  then  used  to  obtain  s(j  -1- 1). 
The  specific  form  of  the  finite  difference  approximation  in  ^{j),  j  €  J[2,  li],  will  be 
discussed  at  a  later  point  in  this  section.  Also,  because  of  the  definition  of  k{j  -h  1), 
after  each  p(j  1)  is  calculated,  it  must  be  checked  to  determine  whether  it  is  such 
that  0.5  <  p{j  + 1)  <  1.5.  If  p(j  +  1)  is  calculated  to  be  such  that  p(j  +  1)  <  0.5, 

1)  must  be  set  to  k{i  + 1)  =  k{j)  -  1  and  p(j  +  1)  must  be  updated  to 

p(j  +  1)  -+  p(j  +  1)  +  1.  If  p(i  +  1)  is  calculated  to  be  such  that  p(>  +  1)  >  1.5, 

k{j  -I- 1)  must  be  set  to  k{j  + 1)  =  k{j)  +  1  and  p(j  +  1)  must  be  updated  to 

P(i  + 1)  p(j  +  1)  -  1. 

At  each  time  step,  once  p(j  + 1),  k(j  + 1),  and  s(j  + 1)  are  determined,  and 
p(j  +  1)  is  updated  if  necessary,  the  displacements  at  the  nodes  are  determined. 
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For  each  j  6  -  1],  the  displacements  u{i,j  +  1),  i  e  +  1)  —  1], 

and  tii(i,>+l),  i  e  I[A:(j  +  1)  +  l,n  —  1],  are  determined  from  centeied-type 
difiference  equations.  It  will  also  be  assumed  in  the  following  that  |s(j  +  1)  -  s(j)l  < 
^/4,  V  j  6  X[0,/i].  This  is  done  so  that  u(i,j  -1),  u(t,j),  and  + 1) 
in  the  difference  equations  representing  Equation  (72)  all  correspond  to  phase  1, 
and  Ui(i,j  —  1),  ui(i,j),  and  ui(i,j  +  1)  in  the  difference  equations  representing 
Equation  (7.3)  all  correspond  to  phase  2.^ 

To  obtain  the  initial  displacements  at  the  nodes,  we  use  the  initial  conditions 
given  by  (7.9)i  and  (7.10)i.  We  assume  that  these  initial  displacements  are  continuous 
V  X  €  (0, 1),  satisfy  the  boundary  conditions  given  by  (7.7)  and  (7.8),  are  such  that 
the  linear  momentum  jump  condition  given  by  (7.5)  is  satisfied,  and  have  a  first  mode 
type  mode  shape.  In  particular,  we  assume  that 

h(x)  =  eoi ,  (8.4) 


for  0  <  X  <  S(0),  and 

fc.(i)  =  --r{(i  -  if  -  (S(0)  -  1)'}  +eoS(0),  (8S) 

for  1(0)  <  X  <  1.  The  values  of  the  displacements  at  the  next  time  increment 
can  be  approximated  by  the  Taylor  series  expansion  in  time  of  the  displacements. 
In  particular,  using  the  initial  conditions  given  by  (7.9)  and  (7.10),  the  first-order 
2q)proximation  of  u  and  ui  at  f  =  T  are 

u(x,  T)  =  h(x)  -i-  g{x)T ,  (8.6) 

*  Note  that  thia  is  coaaMtetu  with  the  aasumpiioo  that  the  mapiinirtet  of  the  leims  in  in  die  linear  momrntmn 
jump  cooditioa  are  negligibie  in  comparison  to  the  mapiitudea  of  the  fint-oider  lenns  in  that  jump  condition. 


18 


for  0  <  X  <  S(0),  and 

ui(x,  T)  =  hi(x)  +  gi{x)T ,  (8.7) 

for  1(0)  <  X  <  1,  respectively.  It  is  assumed  that  the  initial  velocity  distiibutimi 
results  in  a  continuous  displacement  at  the  phase  boundary  at  time  t  =  T,  results 
in  the  boundary  conditions  given  by  (7.7)  and  (7.8)  being  satisfied  at  time  i  = 
results  in  the  linear  momentum  jump  condition  given  by  (7.5)  being  approximately 
satisfied  at  time  i  =  T,  and  has  a  first-mode  type  velocity  profile.  Such  an  initial 
velocity  distribution  is  given  by 


for  0  <  X  <  1(0),  and 


^(x)  =  vqx  , 


(8.8) 


“  2£(7(r~i3<<^  ■ 

for  1(0)  <  X  <  1.  This  initial  velocity  distribution  was  used  in  (8.6)  and  (8.7),  which 
were  then  used  to  obtain  the  displacements  at  the  nodes  at  f  =  T.^ 

For  j  6  —  1],  u(0,y  -1- 1)  is  obtained  firom  the  following  equation  which 

represents  the  fixed  boundary  condition  given  by  (7.7): 

u(0,i-|-l)  =  0.  (8.10) 

V  j  €  -  1).  We  can  obtain  u(*,j  -|- 1),  for  t  6  1(1,  -I- 1)  -  2]  at  each 

j  €  X[l,  h  -  1],  from  a  finite  difference  approximation  of  Equation  (7.2)  which  uses 
centered  difference  equations  for  equally  spaced  nodes  to  approximate  each  term  in 

^  We  aott  that  i(T)  is  not  a  given  consiaat  in  the  proUem  and  Ifaeiefote  cannot  be  naed  if  an  analytical  aohnion 
was  to  be  obtained.  In  this  case,  we  could  use  s(0)  instead  of  i(T)  in  (8.9)  as  an  appnniniatioiL  For  the 
numeiical  method,  however,  we  can  nse  i(T). 
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that  equation.  In  the  following,  centered  difference  equations  for  equally  spaced 
nodes  will  be  referred  to  as  standard  centered  difference  equations.  The  resulting 
difference  equation  that  is  used  to  obtain  u(t,  j  +  1)  is 

«(*>J  +  1)  =  2u(i,  -  l)+a{u(x  +  1,  j)  -  2u(i,  j)  +  u(:  -  1,  j)}  ,  (8.11) 

for  i  €  I[l,  k{j  +  1)  -  2]  at  each  j  €  1(1,  /i  —  1),  where  a  =  {T/hf.  The  difference 
equation  given  by  (8.11)  is  commonly  used  for  the  wave  equation  [2].  It  is  also 
well  known  that  a  numerical  routine  using  this  difference  equation  is  numerically 
stable  if  and  only  if  a  <  1,  and  that  the  numerical  error  increases  as  a  decreases 
from  1  [2].  We  can  obtain  + 1),  for  *  €  I[k(j  +  l)  +  2,n-l]  at  each 

j  €  J[l,  /)  -  1],  from  a  finite  difference  approximation  of  Equation  (7.3)  which  uses 
standard  centered  difference  equations  to  s^proximate  each  term  in  that  equation.  The 
resulting  difference  equation  that  is  used  to  obtain  + 1)  is 

«i(*i  J  +  1)  =  2ui(i,i)  -  u,(t,  J  -  1)  +  Ea{ui{i  +  1,  j)  -  2u,(t,  j)  +  u,(i  -  1,  j)} 

+  7o{s0‘  +  1)  -  2s(j)  +  s{j  -  1)} , 

(8.12) 

for  i  G  I[k{j  +  1)  +  2,n  —  1]  at  each  j  G  Tll,li  —  1].  For  j  G  I[l,  h  —  1].  we  can 
obtain  u(n,  j  + 1)  from  an  O(h^)  difference  equation  representing  the  free  boundary 
condition  given  by  Equation  (7.8).  This  difference  equation  is 

wi(«,  J  +  1)  =  3  {4wi(«  -  1,  J  +  1)  -  -  2,  J  +  1)}  ,  (8.13) 

for  each  j  G  1(1,  /i  -  1]. 

For  the  displacements  u{k{j  +  1)  —  l,i  +  1)  and/or  ui(ib(j  +  1)  +  l,i  + 1), 
finite  difference  methods  using  q{j)  wrill  be  used.  For  the  derivation  of  the  difference 
equations  for  these  displacements,  we  first  note  that  the  distance  separating  s(j) 
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from  its  nearest  nodes  is  not  equal  to  h.  Because  of  this,  difference  equations  for 
unequally  spaced  nodes  must  be  used  to  represent  and/or  These  difference 
equations  will  be  derived  from  the  second-degree  Lagrange  interpolating  polynomials 
for  u  and/or  ui  near  the  interface.  For  a  function  g(x),  its  second-degree  Lagrange 
polynomial,  denoted  by  P(x),  is  given  by 


P{x) 


(x  -  Xi)(x  -  X2) 
(xo  -  Xi)(xo  -  12) 


qM  + 


+ 


(x  -  xo)(x  -  X2) 
(xi  -  xo)(xi  -  X2) 
(x  -  Xo)(x  -  Xi) 
(X2  -  Xo)(X2  -  Xi) 


s(®i) 


(8.14) 


where  xq,  xi,  X2,  <7(xo),  5(xi),  and  ff(x2)  are  given  (see  [2]).  The  first  derivative 
of  P(x)  «  g(x)  is 


dP(x) 

dx 


X  —  Xi  +  X  —  X2 
(xo  -  xi)(xo  -  X2) 


vM  + 


+ 


X  —  Xo  +  X  —  X2 

(xi  —  Xo)(xi  —  X2) 
X  —  Xo  ~h  X  —  Xi 

(12  -  Xo)(X2  -  Xi) 


9(X2), 


(8.15) 


and  the  second  derivative  is 


dx^ 


_ 2 

(xo  —  xi)(xo  —  X2) 


sM  + 


+ 


(xi  -  xo)(xi  -  X2)^^*^^ 

—  ■  —  ^ - o(x2) 

(x2-xo)(^2-ari)^^ 


(8.16) 


If  we  choose  x  =  xo,  xi,  or  X2  in  Equations  (8.15)  and  (8.16),  the  errcvs  in  these 
equations  are  approximately  0{<P),  where  d  is  the  maximum  distance  between  any 
of  these  three  points.^  Also,  we  note  that  when  xo,  xi,  and  X2  ate  equally  spaced, 
(8.16)  reduces  to  a  standard  centered  difference  equation  of  the  form  that  was  used 
in  (8.11)  and  (8.12). 

*  See  [2]  for  a  mote  detailed  discussion  of  the  error  that  is  involved  in  the  Lagnuige  interpolating  polynomial 
and  its  derivatives.  Also,  difference  equations  of  this  type  are  used  in  some  finite  difference  methods  for  Stephan 
problems  (see  [3]). 
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Case  I:  Consider  the  case  where  0.5  <  p(i  +  1)  <  1.5,  before  being  updated  (Hgure 
la).  For  this  case,  k{j  +  1)  =  k(j),  and  for  the  calculation  of  u(k(j  +  1)  —  1,  j ■  +  1), 
a  difference  equation  of  the  form  (8.16)  will  be  used  to  approximate  ati  =  jT 
in  Equation  (7.2).  Using  this  difference  equation  and  a  standard  centered  difference 
equation  for  the  second  time  derivative  term  in  Equation  (7.2),  we  obtain 

«(^0'  +  1)  -  1,  J  +  1)  =  2u(k(j  +  1)  -  l,i)  -  u{k{j  +  1)  -  1,  j  -  1) 


ru(fc(j)-2,j)  u(fc(j)-l,i)  qjj)  \ 

I  1+pO)  pU)  [1 + 

Similarly,  for  the  calculation  of  ui{k{j  + 1)  +  l,i  +  1),  a  difference  equation  of  the 
form  (8.16)  will  be  used  to  approximate  at  t  =  jT  in  Equation  (7.3).  Using  this 
difference  equation  and  standard  centered  difference  equations  for  the  second  time 
derivative  terms  in  Equation  (7.3),  we  obtain 

«i(Mj  +  1)  +  l,i  +  1)  =  2u,(A;(j  +  1)  +  1,  j)  -  ui{k{j  + 1)  + 1,  j  -  1) 


+  2f;o 


{f 


9(j) 


P(j)][3  -  p(i)] 


MHi)  + 1, j)  MKi)  +  2, i) 


2  --  pij) 


3-p(j) 


} 


(8.18) 


+  7o{50'  +  1)  -  2s(j)  +  s(j  -  1)}. 

We  note  that  when  p(j)  =  1,  Equations  (8.17)  and  (8.18)  reduce  to  Equations  (8.11) 
and  (8.12),  respectively. 

Case  U:  Consider  the  case  where  p(J  +  1)  <  0.5,  before  being  updated  (Hgure 
lb).  For  this  case,  k(j  + 1)  =  k(j)  —  1  and  p{j  + 1)  -♦  p(j  + 1)  + 1.  Additionally, 
in  this  case,  we  can  use  the  difference  equation  given  by  (8.11)  for  the  calculatitm 
of  u(k(j  + 1)  -  l,i  +  1).  However,  for  the  calculation  of  Ui(Jb(j  + 1)  +  1,  j ■  + 1),  a 
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difference  equation  of  the  form  (8.16)  will  be  used  to  approximate  and  standard 
centered  difference  equations  will  be  used  to  approximate  the  two  time  derivative 
terms  in  Equation  (7.3).  The  resulting  difference  equation  is 

Ui(Ki  +  1)  +  li i  +  1)  =  2ui(ik(i  +  1)  +  1,  j)  -  Ui(fc(ji  +  1)  +  1,;  -  1) 


+  2Eoc 


{ 


qU) 

[1  -  p(i)]l2  -  p(i)] 


l-p(i)  2-p{j)  j 


(8.19) 


+  7o{a(j  +  1)  -  2s{j)  +  sU  -  1)}. 


Case  III:  The  last  case  that  can  occur  is  p{j  + 1)  >  1.5,  before  being  updated 
(Rgure  Ic).  In  this  case.  k(j  +  1)  =  k{j)  +  1  and  p(j  +  1)  p(j  +  1)  -  1.  In  this 
case,  for  the  calculation  of  u{k{j  4- 1)  —  1,;  + 1),  a  difference  equation  of  the  form 
(8.16)  will  be  used  to  approximate  ^  and  a  standard  centered  difference  equation 
will  be  used  to  approximate  the  second  time  derivative  term  in  Equation  (7.2).  The 
resulting  difference  equation  is 

«(^0'  +  1)  -  1,7  +  1)  =  2uik{j  + 1)  -  1,;)  -  u{k{j  +  1)  -  1,;  -  1) 


fu(fc(j)-l,j)  u(fc(j),j)  qjj)  ] 

1  pO)  pO)-!  p(j)Ip{7)'-i]J  ’ 

(8.20) 


For  this  case,  we  can  use  the  difference  equation  given  by  (8.12)  for  the  calculation 
of  ui(k{j  +  1)  +  1,7  +  1). 


Once  s{j  + 1)  and  the  displacements  at  the  nodes  t  €  J[0,  k{j  +  1)  —  1]  and  i  e 
J[ib(j  +  1)  +  1,  n]  have  been  determined,  we  determine  the  displacement  q(j  + 1)  at 
the  interface  from  the  difference  equation  representing  the  linear  momentum  jump 
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condition.  Using  difference  equations  of  the  form  (8.15)  to  approximate  the  spatial 
derivatives  in  Equation  (7.5),  we  obtain 


sO'  + 1) 


IP- 


E^p  —  5]  2p  1 

-pi(3-p]“(r+ 


— I  (- 

-pIpj  (I 


+  1) 


-  i,i  + 1)  -  fif— + i,i  + 1) 

P  2  -p 


(8.21) 


-E^^^ui{k  +  2,j  +  1)| , 

where  k  =  k{j  +  1)  and  p  =  p( j  + 1). 

The  last  displacement  that  must  be  calculated  at  each  time  increment  is  the 
displacement  at  the  node  k(j  +  1).  This  displacement  will  be  calculated  using 
a  second-degree  Lagrange  interpolating  polynomiaL  In  particular,  if  k(j  +  l)h  < 
s{j  4- 1),  we  calculate  u(Jt(j  4-  l),i  4- 1)  from 


“(t-J  + 1)  =  -  2,  j  + 1)  +  - 1,  j  +  1)  +  ■  <*■“> 

1  4-  p  p  11  4-  pjp 


where  k  =  k{j  4- 1)  and  p  =  p(j  4- 1),  and  if  k{j  4-  l)h  >  s{j  4- 1),  we  calculate 
4-  l),j  4- 1)  from 


(»3S) 


where  k  =  k{j  4- 1)  and  p  =  p{i  4- 1). 


The  last  quantity  that  is  calculated  at  each  tinae  step  is  i{j  4-  !)•  Recall  that 
this  quantity  is  used  in  the  calculation  of  p{j  -f- 1)  at  the  next  time  step.  For  the 
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calculation  of  $(j  +  1),  j  €  -  1],  a  finite  difference  equation  of  the  form 

(8.1S)  is  used  for  the  approximation  of  The  resulting  finite  difference 

equation  for  +  1)  is 

«(i  + 1)  =  -  -  2.y  + 1)  -  -  l.i  + 1) 


2p-n 
[1  +  p]p 


9(j  + 


(8.24) 


for  j  €  T[l,/i  -  1],  where  k  =  k(j  +  1)  and  p  =  p(j  +  1). 

The  numerical  routine  that  was  discussed  above  allows  for  the  phase  boundary  to 
pass  over  nodes  other  than  the  node  k{j).  However,  for  a  problem  that  uses  a  kinetic 
relation  of  the  form  given  by  (7.6),  most  values  of  the  material  coefficimts  that  are 
consistent  with  the  assumption  that  the  second-order  terms  in  the  linear  momentum 
jump  condition  are  negligible  in  comparison  to  the  first-order  terms  in  that  equation 
and  most  initial  conditions  that  produce  infinitesimal  initial  stains  will  result  in  the 
phase  boundary  staying  within  the  interval  between  the  nodes  k{j)  -  1  and  k{j)  + 1 
for  all  f  €  r.  In  this  case,  a  simplified  numerical  routine  can  be  used  where  k(j)  has 
the  same  value  at  each  time  increment,  p{j  +  1)  never  needs  to  be  updated  (in  the 
sense  that  it  was  updated  in  Clases  II  and  III),  and  only  Case  I  for  the  calculation  of 
the  displacements  near  the  interface  needs  to  be  considered. 


9.  The  Free  Vibrations  and  Damping  Properties 

In  this  section,  the  free  vibrations  of  the  two-phase  bar  that  were  determined 
from  the  finite  difference  method  that  was  discussed  in  the  previous  section  are 
discussed.  As  expected,  the  response  of  the  two-phase  bar  to  the  initial  conditions 
given  by  (8.4),  (8J),  (8.8),  (8.9),  and  (7.11)  has  a  decaying  oscillatory  form.  In 
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particular,  the  position  of  the  phase  boundary  oscillates  as  time  progresses  and  decays 
to  a  new  position  that  has  a  distance  and  direction  from  its  initial  position  that  is 
proportional  to  the  magnitude  and  “direction”  of  its  initial  conditions  (Figures  2-4). 
The  displacements  also  Irave  a  decaying  oscillatory  form,  and  they  go  to  zero  as  time 
goes  to  infinity  (Figures  5-9).  The  mode  shape  of  the  bar  during  these  free  vibrations 
has  a  first  mode  type  form  (Figure  10)^  This  is  most  likely  a  result  of  the  fact  that 
first  mode  type  initial  conditions  were  given. 

9.1.  The  Damping  Behavior 

The  damping  of  the  bar  was  studied  as  and  70  were  varied.  In  Figure  11, 
a  plot  of  the  settling  time  versus  v  is  presented. Here,  the  settling  time  is  defined  as 
the  nondimensional  time  necessary  for  the  amplitude  of  ui{L,t)  to  become  less  than 
10~^.  As  the  settling  time  decreases,  it  is  said  that  the  damping  of  the  bar  increases. 
As  expected,  as  i>  decreases,  the  damping  of  the  bar  increases.  This  is  primarily  a 
result  of  the  fact  that  for  a  given  amount  of  strain  at  the  interface,  the  nominal  phase 
boundary  velocity  increases  as  u  decreases.  Consequently,  as  i>  decreases,  there  is 
more  motion  of  the  phase  boundary  in  a  given  interval  of  time,  which  results  in 
more  energy  being  dissipated  in  that  interval  of  time.  This  increase  in  damping  as 
v  decreases  is  also  displayed  in  Figures  2-7.  From  Figure  11,  it  appears  that  the 
“frequency”  of  oscillation  does  not  significantly  depend  on  u.  In  particular,  it  can  be 
observed  from  this  figure  that  as  i/  is  varied  the  settling  time  remains  constant  over 
an  interval  of  i/,  and  when  the  settling  time  changes,  it  does  so  discontinuously.  This 

*  In  Figure  10, «  vs  t  is  plotted  for  0  <  x  <  i{i)  ind  i,  vs  tit  plotted  for  i(t)  <  z  <  1.  The  shape  defonnanon 
for  phase  2  is  not  plotted. 

In  this  hgure,  the  settling  tiine  versus  v  is  not  presented  for  0  <  v  <  0.09  because  the  niuuetical  routine  is 
unstable  for  these  values  of  v  when  E  =  1.15,  -n  =  0.1,  T  =  0.01,  and  h  =  0.02.  This  will  be  dtmiMed 
further  in  the  next  section.  Also,  in  this  6gure  and  in  Figures  12  and  13,  the  initial  conditions  that  were  used  are 
eg  =  0.001,  eg  0.001,  and  a«  s:  0.5. 
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is  a  reflection  of  the  fact  that  as  i>  decreases  the  amplitude  of  oscillation  decreases, 
but  the  “frequency”  of  oscillation  does  not. 

In  Figure  12,  the  settling  time  versus  E  is  presented.  From  this  figure  it  can  be 
concluded  that  as  E  goes  to  zero,  the  damping  goes  to  zero  (i.e.  the  settling  time  goes 
to  infinity),  and  as  E  increases,  the  damping  increases.  Additionally,  this  increase 
in  damping  levels  off  at  a  relatively  high  level  of  damping  (i.e.  at  a  relatively  short 
settling  time)  after  E  fa  0.15. 

For  the  plots  of  the  settling  time  versus  the  transformation  strain  70,  we  do 
not  vary  70  and  keep  E  and  i>  constant.  This  is  because  E  and  i>  can  remain 
constant  as  70  is  varied  only  if  Ei/E  and  also  change  values.  We  instead  let 

E'  =  Ei/E  and  i/'  =  and  substitute  E  =  E'/{1  +  70)  and  i>  =  v'/'yo  into  the 

difference  equations  of  the  finite  difference  program.  We  then  plot  the  settling  time 
versus  70  while  keeping  E'  and  P'  constant  An  example  of  such  a  plot  is  shown  in 
Figure  13.  From  this  figure  we  can  conclude  that  as  70  goes  to  zero  the  damping  goes 
to  zero,  and  as  the  magnitude  of  70  increases  the  damping  increases.  We  can  also 
observe  from  this  figure  that  as  the  magnitude  of  70  becomes  greater  than  |7o|  »  0.06 
the  damping  becomes  appreciable,  and  for  I70I  less  than  this  value  the  damping  is 
relatively  small.  It  is  interesting  to  note  that  if  one  were  to  define  an  infinitesimal 
transformation  strain  as  one  that  produces  a  vibration  response  that  is  qualitatively  like 
that  produced  as  70  — »  0,  and  if  one  were  to  define  a  finite  transformation  strain  as  one 
that  produces  a  vibration  response  that  is  qualitatively  like  that  produce  by  a  70  that  is 
close  to  or  less  than  -0.5  or  is  close  to  or  greater  than  1,  for  the  values  of  E'  and  v* 
considered,  a  transformation  strain  with  a  magnitude  that  is  less  than  {70 1  »  0.06 
would  be  considered  infinitesimal,  and  a  transformation  strain  with  a  magnitude  that 
is  greater  than  I70I  0.06  would  be  considered  finite.  These  “transition”  values 
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7o  ^  -0.06  and  70  ^  0.06  separating  the  infinitesimal  transformation  strains  from 
the  finite  transformation  .strains  are  probably  much  smaller  than  most  might  have 
guessed  beforehand.  This  also  underlines  the  importance  of  treating  a  transformation 
strain  that  is  not  truly  infinitesimal  as  a  finite  strain,  for  at  least  vibration  problems. 

9,2.  Instabilities  of  the  Numerical  Routine 

As  mentioned  previously,  there  are  some  instability  problems  with  the  numerical 
routine  for  values  of  v,  £,  and  70  outside  of  a  certain  region  of  the  parameter  space. 
In  some  cases,  these  instabilities  resemble  those  of  highly  damped  systems,  which 
are  sometimes  referred  to  as  stiff  systems  in  the  numerical  methods  literature  [2]. 
In  the  remaining  cases,  the  loss  of  stability  of  the  numerical  routine  resembles  that 
of  a  standard  centered  difference  numerical  routine  for  a  wave  equation  when  the 
coefficient  multiplying  the  term  representing  the  spatial  derivative  becomes  greater 
than  one.  In  both  cases,  however,  the  values  of  u,  £,  and  70  where  the  numerical 
routine  loses  stability  depends  on  the  values  of  T  and  h  that  are  used.  For  example, 
for  h  =  0.02,  E  =  1.15  and  70  =  0.1,  the  numerical  routine  loses  stability  as  i> 
is  decreased  at  i/  w  0.22  when  T  =  0.018,  and  at  P  w  0.09  when  T  =  0.01.  The 
loss  of  stability  in  both  of  these  cases  resembles  that  of  stiff  systems.  In  fact,  for 
the  latter  case,  the  displacements  reach  their  settling  time  in  almost  one  half  of  one 
“cycle”  of  oscillation.  One  should  note  however  that  as  i>  gets  close  to  zero,  the 
assumption  that  1/  is  such  that  the  magnitude  of  is  negligible  in  comparison  to 
the  magnitudes  of  the  first-order  terms  in  Equation  (4.7)  becomes  less  valid.  For 
h  =  0.02,  u  =  0.5,  and  70  =  0.1,  the  value  of  E  beyond  which  the  numerical 
routine  is  unstable  is  E  as  1.23  when  T  =  0.018,  and  E  »  4.01  when  T  =  0.01. 
Both  of  these  values  of  E  correspond  to  E{T/hf  »  1  (recall  that  this  term  appears 
in  the  finite  difference  approximation  of  the  forced  wave  equation  for  phase  2  given 
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by  (8.12)).  For  T  =  0.018,  h  =  0.02,  E'  =  1.265,  and  u'  =  0.05,  the  numerical 
routing  is  stable  for  0.021  <  70  <  0.152.  The  loss  of  stability  at  70  «  0.021 
corresponds  to  E{T/hf  w  1,  and  the  loss  of  stability  at  70  «  0.152  resembles  that  of 
a  stiff  system.  For  T  =  0.01,  h  =  0.02,  E'  =  1.265,  and  v'  =  0.05,  the  numerical 
routing  is  stable  for  -0.232  <  70  <  0.256.  For  this  case,  the  loss  of  stability  at  both 
7o  w  0.256  and  70  w  -0.232  resembles  that  of  stiff  systems. 
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